Effect of stress-triaxiality on void growth in dynamic fracture of metals: a molecular 

dynamics study 
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The effect of stress-triaxiality on growth of a void in a three dimensional single-crystal face- 
centered-cubic (FCC) lattice has been studied. Molecular dynamics (MD) simulations using an 
embedded-atom (EAM) potential for copper have been performed at room temperature and using 
strain controlling with high strain rates ranging from 10 7 /sec to 10 10 /sec. Strain-rates of these 
magnitudes can be studied experimentally, e.g. using shock waves induced by laser ablation. Void 
growth has been simulated in three different conditions, namely uniaxial, biaxial, and triaxial expan- 
sion. The response of the system in the three cases have been compared in terms of the void growth 
rate, the detailed void shape evolution, and the stress-strain behavior including the development 
of plastic strain. Also macroscopic observables as plastic work and porosity have been computed 
from the atomistic level. The stress thresholds for void growth are found to be comparable with 
spall strength values determined by dynamic fracture experiments. The conventional macroscopic 
assumption that the mean plastic strain results from the growth of the void is validated. The evolu- 
tion of the system in the uniaxial case is found to exhibit four different regimes: elastic expansion; 
plastic yielding, when the mean stress is nearly constant, but the stress-triaxiality increases rapidly 
together with exponential growth of the void; saturation of the stress-triaxiality; and finally the 
failure. 

PACS numbers: 61.72.Qq, 62.20.Mk, 62.20.Fe, 62.50.+P 
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I. INTRODUCTION 



Ductile fracture of metals commonly occurs through 
the nucleation, growth and coalescence of microscopic 
voids. 1 Much can be learned about the ductile fracture 
through the study of these voids. A particularly interest- 
ing case is the dynamic fracture of ductile metals^^^ 
in which the strain rates are so high that processes such 
as diffusion operating on relatively long time scales may 
be neglected, while inertial effects become relatively im- 
portant. Void growth is driven by the need to relax ten- 
sile stress that builds up in the system, and to mini- 
mize the associated elastic energy. The material around 
a void deforms plastically in order to accommodate the 
void growth. Naturally, the plastic deformation results 
from a local shear stress, which may arise from the ap- 
plied stress, but it also may arise from the stress field 
of the void even if the applied stress is hydrostatic. So 
the expectation is that the evolution of the plastic zone, 
and hence the growth of the void, is influenced by the 
degree of stress triaxiality; i.e., the ratio of the mean 
(hydrostatic) stress to the shear stress. It is this rela- 
tionship that we study here by varying the triaxiality 
of the loading. In particular, we conduct simulations 
in which one, two or three directions of the system are 
expanded, producing a state of uniaxial, biaxial or tri- 
axial strain, respectively. Variation in the triaxiality of 
the strain causes variation in the triaxiality of the stress 
state, where it should be noted that uniaxial (biaxial) 
strain does not imply pure uniaxial (biaxial) stress. 

Besides dynamic crack propagation experiments, dy- 
namic fracture can be measured for instance in shock 
physics or spallation experiments, to which the simula- 



tions performed here are compared. Various techniques 
are employed to generate the shock waves: Hopkinson 
bar, gas gun, high explosives, and laser ablation. With 
Hopkinson bar the strain-rates i usually are of the or- 
der I0 2 — 10 4 /sec, in gas gun of 10 5 /sec, with high- 
explosives even higher strain-rates can be produced, and 
with lasers strain-rates exceeding 10 7 /sec are attained. 
In a gas gun for instance the fracture results from es- 
sentially one-dimensional shock loading. Two compres- 
sive shock waves are generated by the impact of a flier 
on a metal target, propagate away from each other, re- 
flect from opposite free surfaces becoming tensile release 
waves and finally come into coincidence again. If the 
combined tensile stress exceeds the rupture strength of 
the material, the metal fails, after some incubation time, 
producing a fracture surface. In strong shocks, a scab of 
material may spall from the back side of the target and 
fly off. Spallation experiments^ for single and polycrystal 
copper report spall strength values of a* ~ 3 — 4 GPa at 
strain-rates e ~ 2 — 3 x I0 5 /sec, and scaling between the 
spall strength and strain-rate, a* ~ £ - 2 . 

In this study of dynamic fracture in ductile metals at 
high strain-rates (10 7 - 10 10 /sec) we have concentrated 
on void growth starting from a single crystal copper lat- 
tice containing an infinitely weakly bound inclusion or a 
pre-existing nanoscale void. The lattice is initially free 
of other defects. We have focused on the effect of stress- 
triaxiality on void growth. In some fracture experiments, 
for example in necking and cup-cone fracture^ the uniax- 
ial strain produces a stress state that transitions rapidly 
to triaxial state due to the plastic flow during the course 
of loading. It is during the triaxial phase that void growth 
and failure take place. Because of the connection with 
shock experiments, the stress triaxiality study done here 
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is carried out using strain control, and it is the strain 
that is varied from uniaxial to biaxial to triaxial. 

Much of the damage modeling of metals has been 
carried out at mean-field or continuum level based on 
constitutive theories. The continuum models concen- 
trate especially on two areas: macroscopic crack growth 
phenomenon^ and studies of porosity, i.e. behavior 
of an array of voids, at sub-grain level during load- 
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area for example the locus of yield surfaces in stress 
space has been studied, which is related to the ques- 
tion of the effect of stress-triaxiality studied here. Of 
the void growth studies especially the Gurson modeli^ 
is commonly used to model cavitation (the development 
of porosity) at the sub-grid level in what is termed as 
damage modeling. These continuum calculations often 
assume that the matrix material, where the voids are 
embedded, is elastically rigid and plastically incompress- 
ible, and the dilation of the void-matrix aggregate is com- 
pletely due to the void growth. Of particular interest, 
and relevance in terms of this study, is a single crys- 
tal plasticity study of void growthiSl The calculations 
are typically done by determining approximate solutions 
for integrals of incremental equations of virtual work us- 
ing the finite element method. Continuum modeling has 
been used to study some of the phenomena addressed in 
this paper such as the effect of triaxiality on void growth 
and void shape changes^ The validity of the approxi- 
mate solutions of incrementals limits the strain-rates to 
be rather low compared to the strain-rates used in this 
study. 

In order to characterize the void growth not only with 
macroscopic quantities and at the continuum level, but to 
investigate what happens at the atomistic level, we have 
employed molecular dynamics (MD) simulations. MD 
simulations enable us to see what the effects are on the 
void surface at the single atom level, when it grows while 
the total system yields. This Article presents work, which 
is an extension to the work done earlier by some of this 
Article's authors of void growth in a single crystal copper 
with hydrostatic loading, void nucleation and growth in 
single and polycrystalline CO pperi22i2Si2ii224^iM Molecu- 
lar dynamics simulations of void growth in single crystal 
copper have also been conducted by other groups for slab 
geometries of interest in the semiconductor chip metaliza- 
tion problemii&Si The effectively two-dimensional, thin 
film systems are in contrast to three dimensional bulk 
systems studied here. In some cases MD simulations of 
plasticity and crack propagation have been as large as 
billion atoms 4£ 

This Article is organized as follows. It starts in Sec- 
tion[n]with an overview of the MD method used and the 
simulations which have been carried out in this study. 
Exploration of the results of the simulations starts in Sec- 
tion IIIII by the study of the mean or hydrostatic stress 
versus strain as well as the deviatoric part of the stress 
tensor, von Mises stress, which is used to measure the 
shear stress, and stress-triaxiality versus strain for all 



the simulated strain-rates and modes of expansion. Sec- 
tion [^3 concentrates on the macroscopic plastic quanti- 
ties, such as mean and equivalent plastic strain, plastic 
work and its relation to the temperature. The evolution 
of the void in terms of its growth and shape changes is 
studied in Section Section IVll summarizes the results 
and compares different measured quantities with each 
other concentrating on one of the simulations, uniaxial 
strain with strain-rate 10 8 /sec. The Article is concluded 
with discussions of the results and suggestions for future 
studies in Section IVIll 



II. METHOD AND SIMULATIONS 

A. Strain-controlled Molecular Dynamics 

In this atomistic-level study of void growth, the sim- 
ulations have been done using empirical embedded-atom 
(EAM) potentials in classical molecular dynamics^ fol- 
lowing the scheme developed earlier ,22i2fli2i The copper 
EAM potential we have used is due to Oh and John- 

SOQpM 

The system, in which the simulations are done, 
is a three-dimensional single-crystal face-centered-cubic 
(FCC) lattice in a cubic box with {100} faces. Periodic 
boundary conditions are used in all the three directions 
so that there are no free boundaries in the system apart 
from the void. Equivalently the system can be imagined 
to consist of an infinite periodic array of voids. Note 
that periodic boundary conditions have also been used 
in continuum models of void growth, but in the contin- 
uum modeling of void growth in isotropic materials the 
calculations are done in a reduced cell, which exhibits 
one quarter of the box in two dimensions, and one eighth 
in three dimensions, and the behavior of other areas are 
derived from the symmetries. We use the full cubic box 
because the cubic symmetry present in the continuum is 
broken in MD at finite temperature, and processes such 
as dislocation nucleation at the void surface would be 
over-constrained in a reduced simulation box. 

In the simulations, the system is brought to thermal 
equilibrium at room temperature, T — 300 K, with a 
commonly used thermostalssi and at ambient pressure, 
P ~ MPa, keeping the volume constant. After that a 
spherical void is cut in the middle of the system, later 
the thermostat is turned off, and the dilational strain 
is applied uniformly with a constant strain-rate i. The 
removal of the atoms in the spherical region may be con- 
sidered to simulate the instantaneous separation of the 
matrix material from an infinitely weakly bound inclu- 
sion. The uniform expansion in these strain-controlled 
simulations is applied through rescaling the coordinates 
as in the Parrinello-Rahman method.— Technically the 
three Cartesian coordinates of the atoms are rescaled to 
the unit-box, each coordinate S a £ [0,1). When calcu- 
lating the forces and velocities, as well as updating the 
new positions of the atoms, the unit-box is multiplied by 
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a diagonal scaling matrix 7i — {l x ,l y ,l z }, where I's are 
the simulation box's side lengths, to compute the true 
positions of the atoms, 

x = HS. (1) 

This scaling matrix TL is updated each time-step, when 
the load is applied, by multiplying the initial matrix TLq 
with the sum of the unit matrix and the strain matrix 
£ = t£, 

TL(t) =TL (l + t£). (2) 

For our purposes the strain-rate matrix £ is always diag- 
onal, since neither rotation nor simple shear type strains 
are studied. In the triaxial case all the terms in the diag- 
onal are equal; in the uniaxial there is a single non-zero 
term; and in the biaxial case two of the three diagonal 
terms differ from zero and are equal. Prior to expansion 
the system is cubic, its scaling matrix TLq is diagonal, and 
all the terms are equal and correspond to the equilibrium 
size at ambient pressure. Hence the scaling matrix TL re- 
mains diagonal throughout the simulation, and the strain 
in each case is in a (100) direction. 

In fracture and plasticity simulations the first quan- 
tity to consider is the stress-strain behavior. With the 
strain as an input parameter, here we have to measure 
the stress. In this study of the stress-triaxiality we are 
interested in both mean and shear stresses. Therefore the 
whole stress tensor a a p is needed. The stress tensor (the 
negative of the pressure) can be calculated atomistically 
on each time-step using the virial formula^ 

\ i i j>i ) 

The first term in the stress tensor is the kinetic contri- 
bution of atoms denoted with i and having masses rrii 
and momenta pt. The second term, a microscopic virial 
potential stress, consists of sums of interatomic forces 
of atom pairs (ij) with corresponding distances r^ . It 
should be noted, that here and in the rest of the Article 
i and j denote the atoms, and a and (3 the Cartesian 
coordinates. Note that the thermal stress is included, al- 
though in practice in these simulations it contributes less 
than 1 GPa, less than 10% of the yield stress value, and 
never dominates the changes in stress. 

B. Simulations Performed 



The relatively inexpensive potential used enables us 
to do extensive simulations in time. A single time-step 
takes typically about 40 sec of CPU-time in a system 
with 864 000 atoms in a Linux workstation with Intel 
Xeon 1700MHz processor. The longest calculation re- 
quired 835 050 time-steps corresponding to 5.6 nanosec- 
onds. The time-step was 6.7 femtoseconds. 

As mentioned earlier, in order to study the effects of 
the stress-triaxiality and different modes of expansion on 
the void-growth, we have applied three different types of 
expansion, namely uniaxial, biaxial, and triaxial. The 
strain-rates used for each of the three modes of expan- 
sion are e = 10 10 /sec, 10 9 /sec, 5 x 10 8 /sec, 10 8 /sec, and 
10 7 /sec. For the lowest strain-rate, the MD code was 
parallelized in order to take advantage of massively par- 
allel computers. The parallelization was done using a 
spatial domain-decomposition, and was shown to scale 
nearly linearly up to 128 processors. The parallel code 
was used in the case with 835 050 time-steps mentioned 
above, for example. 

For comparison in the elastic regime, we have also per- 
formed simulations without a void in all three modes 
of expansion. These simulations have been used to de- 
termine the bulk, elastic stress-strain response of the 
EAM copper and hence the elastic constants. With- 
out a void, the system is not so strain-rate and sys- 
tem size dependent, at least up to the point of failure, 
so the so-called "no void" simulations have been per- 
formed with a smaller system size, 45 FCC-cells in each 
direction (364 500 atoms) and at the single strain-rate 
i = 10 9 /sec. A uniaxial study of the 60 3 system size, but 
with a smaller initial void radius of 1.1 nm, was carried 
out with the strain-rate e — 10 8 /sec in order to study the 
void-size dependence. In this case the system with the 
void contains 863 543 atoms. 

It should be mentioned, too, that all of the intermedi- 
ate strain-rate simulations (e = 10 8 /sec and 5 x 10 8 /sec) 
expansion were not started from equilibrium conditions 
at P = MPa, but from systems expanded previously at 
the strain-rate e = 10 9 /sec. These simulations have been 
restarted well before yielding, when the system's behav- 
ior is rate independent, and relaxed for 2000 time-steps, 
or 13.4 picoseconds, without expansion before continu- 
ing the expansion at the intermediate strain rates. The 
energy is conserved during the relaxation in MD simula- 
tions. These restarts have been accomplished at strain 
values e = 4.12 %, e = 2.06 %, and e = 1.72 % in uniax- 
ial, biaxial, and triaxial cases, respectively. 



Typically in the simulations carried out here, the cube 
consists of 60 FCC unit cells in each direction, giving 
864 000 atoms. The equilibrium side-length of such a 
copper system is I = 21.6 nm at room temperature and 
ambient pressure. The radius of the spherical void cut 
from the system, unless otherwise noted, is 0.1 of the 
side-length of the box; thus 2.2 nm. After the void is 
cut, there are 860 396 atoms in the system. 



III. STRESS-STRAIN BEHAVIOR AND 
STRESS-TRIAXIALITY 

Let us begin to explore the results of the MD simu- 
lations by looking at the stress-strain curves. Figure 
shows these curves for each of the modes of expansion at 
all the strain-rates computed. The data from "no void" 
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cases are also plotted. The mean or hydrostatic stress, 



= g Tr a a p, 



(4) 



is plotted to indicate the principal impetus for void 
growth. Note that the strain is the engineering strain, 
defined as the expanded system size divided by the orig- 
inal system size minus one. In the uniaxial and biaxial 
cases the strains are the principal strain values e in the 
direction of the strain, such as e x = e, s y = e z = in the 
uniaxial case, and e y = e z = s, e x = in the biaxial case. 
Hence the strain is the value of a non-zero diagonal term 
of t£ in Eq. (J2J). The mean strains e m are 1/3 and 2/3 
of the plotted uniaxial and biaxial strains, respectively. 
Thus the total volumetric strain-rates are not the same 
in the different expansion modes. In the triaxial cases 
the plotted and the mean strain values are the same. 
The shape of the stress-strain curves do not differ much 
depending on the modes of expansion in these cases, at 
least when plotted versus mean strain. Independent of 
the strain-rate, and whether with or without a void, the 
stress-strain curves lie essentially on top of each other 
during elastic expansion, i.e., the initial smooth behavior 
when the system is still recoverable and has not deformed 
plastically. 

The stress-strain curve starts to deviate from the trend 
of the elastic behavior at a specific, "critical" point which 
we call here a yield point. In other quantities we mea- 
sure a change in behavior happens at a specific point, too, 
and as we shall see later the critical or yield points mostly 
coincide with each other, i.e., their strain values are ap- 
proximately the same independent of from which quan- 
tity we derived it. The same point is also the one when 
the void starts to grow, which is the primary mechanism 
for plasticity in this study. Here we define numerically 
the yield point of the cases with void by comparing their 
stress-strain curves with the reference "no void" curve, 
which behaves elastically beyond the yield points of the 
other cases [cf. Fig. Hfa) inset] . Ultimately the no void 
case does fail by homogeneous nucleation of voids, and 
this is the reason for the drop in the mean stress. There 
is a small offset between cases with a void compared to 
the case without a void due to the elastic relaxation of 
the void. The value of the stress at the yield point in the 
cases with void is lower with lower strain-rate, and thus 
the strain to yield is also lower. In each of the modes 
of expansion, the stress at the yield point for the strain- 
rate e = 10 7 /sec is close to the value to which the higher 
strain-rates converge. Of course, at much lower strain 
rates the physics changes, new mechanisms become ac- 
tive, so this value need not hold for arbitrarily low strain 
rates. However, it is noteworthy that the stress at the 
yield point is not scaling with strain, contrary to the 
experimental finding for the spall strength explained in 
Section [J Overshooting, the phenomenon that the max- 
imum stress is much higher than the stress at the yield 
point, is evident here for the higher strain rates. The 
scaling of the spall strength versus the strain-rate 6 with 



an exponent 0.2 is reproduced here when one compares 
the maximum stress values instead of the stresses at the 
yield point for strain-rates i = 5 x 10 8 /sec, 10 9 /sec, and 
10 10 /sec, since then the exponent is 0.14 — 0.18, lowest 
for the uniaxial case and highest for the triaxial case. 
On the other hand the stress value at the yield point, 
which is at the same time the maximum stress, when 
e = 10 7 /sec is very close to the value of 6-8 GPa the 
spall strength scaling predicted from the lower strain- 
rates mentioned in Section Q] It should be noted also, 
that since we are limited to finite, fairly small, system 
sizes, at late stages of the stress-strain curves, at the 
failure, the data is not realistic anymore. The reason 
is that at the plastic part of the stress-strain behavior 
when the void grows, it also emits dislocations, and in a 
finite system with periodic boundaries, when the disloca- 
tions have traveled long enough, they propagate through 
the boundaries and reenter from the other side. In the 
picture where we have a periodic array of voids in an 
infinite system this means that the voids are so close to 
each other that they start to interact. In reality voids are 
never arranged in a perfect cubic lattice structure and in 
symmetric positions with respect to each other, and thus 
the interactions of the voids in the simulations with their 
periodic images are just an unphysical finite size effect. 

In the shear stress or more precisely in the deviatoric 
part of the stress tensor a e plotted in Fig. a much 
bigger difference is seen between the modes of expansion 
than in the mean stress. For the deviatoric part of the 
stress o~ e we use von Mises stress: 



o- e = [3 J 2 ] 



1/2 



(5) 



where J 2 = ^Tr a' 2 is the second invariant of the stress 
deviator a' a0 = a a p — cr m l42 Thus von Mises stress reads: 
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While the mean stress at the yield point gets a value 
of about a = 5.6 — 6.4 GPa when loaded with strain- 
rate e = 10 7 /sec in each of the three modes of expan- 
sion, von Mises stress has a value of er e = 2.0 GPa and 
0.7 GPa in the uniaxial and biaxial cases, respectively. In 
the triaxial case it should be zero by symmetry, and the 
difference from zero, representing symmetry-breaking ef- 
fects, is small. Thus the loading differences between the 
modes of expansion are quantified in von Mises stress. 
After the onset of plasticity or the void growth, von Mises 
stress gets a value of about <7 e = 0.4 GPa, 0.2 GPa, and 
0.1 GPa, in the uniaxial, biaxial, and triaxial cases, re- 
spectively, independent of the strain-rate, but with sig- 
nificant fluctuations in this regime. In this period also 
dislocations move under the action of the shear stress 
until the stress has dropped to the point that it is no 
longer sufficient to move a dislocation through the for- 
est of dislocations. The final value of the shear stress 



5 



corresponds to the flow stress, and they are close to the 
tensile strength values of copper, 200-400 MPa, quoted 
in the literature^ 

Although von Mises stress captures differences between 
the loading modes quite well, an even better quantity to 
study is the ratio between hydrostatic and shear stresses, 
the stress-triaxiality 



X = oW<7 e 



(7) 



which has been plotted in Fig. |3 In the uniaxial case 
the stress-triaxiality starts from the value x — 3-0 and 
slowly decreases linearly to a value x — 2-8 until the on- 
set of rapid growth at the yield point. After the rapid 
increase the stress-triaxiality saturates at \ — 11.0—16.0. 
The stress-triaxiality in the biaxial case starts with a 
much larger value than in the uniaxial case. It begins 
at x — 6.0 and increases linearly to a value x — 8.0 
at the yield point where it grows rapidly to a value of 
X — 15.0 — 30.0. We have noted the correspondence of 
the final von Mises stress and the flow stress above. Simi- 
larly the spall strength provides an experimental measure 
of the mean stress that can be supported in void growth. 
Thus the stress-triaxiality values can be compared with 
the ratio of the spall and the tensile strength of copper. 
Previously quoted literature values for them are 6-8 GPa 
and 200-400 MPa, respectively, giving for their ratio val- 
ues between 15 and 40, and thus comparable with the 
stress-triaxiality values here. The comparison is not fully 
rigorous, but it provides an indication of how reasonable 
the final stress triaxiality values are in terms of experi- 
ment. Since the stress-triaxiality is the mean stress di- 
vided by von Mises stress, which is equal to zero in the 
triaxial case until the yield point and very small even af- 
ter that, the stress-triaxiality is diverging and therefore 
not plotted here in that case. The stress-triaxiality val- 
ues in the uniaxial and biaxial cases at the elastic part 
of the simulation is compared here also with the values 
one get from the elasticity theory^ 



= 1 Cn - 
X 3 C n 
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S = 1 in the uniaxial case and S = 2 in the biaxial 
case. The literature values for the elastic constants of 
copper are Cn = 168 GPa and C X2 = 121 GPa. 46 Thus 
X = 2.9 and x = 5.8 in the uniaxial and biaxial cases, re- 
spectively, which compare quite well with the simulations 
presented here. When the elastic constants are derived 
from the stress-strain curves as e — > 0, they are close to 
the actual experimental values: Cn ~ 162 GPa, C12 — 
121 GPa and Cn ~ 168 GPa, C 12 ^ 124 GPa in the 
uniaxial and biaxial cases, respectively. 

The critical mean stress, von Mises stress, and the 
stress-triaxiality values where their behaviors start to 
deviate compared to the elastic ones, or what we call 
yield points, are summarized in Tabic [i] for the strain- 
rates e = 10 9 /sec, 5 x 10 8 /sec, 10 s , and 10 7 /sec of the 
principal strains. In the case of the highest strain-rate 



e = 10 10 /sec, the shapes of the stress-strain curves are 
so much rounded due to over-shooting, that there is no 
clear point, where the stress-strain curve deviates from 
the elastic behavior, and thus our definition of the yield 
point is no longer suitable. In comparing the mean strain 
values e m (as in the Table) at the onset of plasticity for 
a particular strain rate, one finds that the uniaxial ex- 
pansion always starts to yield at the least strain, and the 
hydrostatic expansion, at the greatest strain. There are 
two effects that contribute to the increase in the plas- 
tic threshold as the triaxiality increases. First, the shear 
component of the applied stress contributes to the re- 
solved shear stress and lowers the threshold for hetero- 
geneous nucleation of dislocations at the void surfaced 
And second, the volumetric strain-rate is lowest in the 
uniaxial case and the highest in the hydrostatic case. 
Strain-rate hardening then leads to an increase in the 
stress value at the onset of plasticity as the triaxiality in- 
creases. The difference between the critical strain values 
when defined as when a behavior deviates from the elas- 
tic behavior is nearly negligible and thus independent 
of whether one uses the criterion from the hydrostatic 
stress, von Mises stress or stress-triaxiality curves. The 
differences reflect mainly the difficulties in defining the 
point what we call the yield point. However, we will see 
later that if the mean stress and von Mises stress start to 
deviate from the elastic behavior with the same ratio as 
they have during elastic expansion, the stress-triaxiality 
may deviate a bit later than the other quantities. We 
shall see later, too, that the onset of plasticity defined 
from these quantities is very close to where the void starts 
to grow. 



IV. PLASTIC STRAIN AND PLASTIC WORK 

After sufficient expansion, the system yields and the 
mean stress is observed to drop with respect to the elas- 
tic response. Then as the simulation box continues to 
expand, the stress remains roughly constant until the 
precipitous drop at final failure. In the region of in- 
creasing strain but roughly constant mean stress, most of 
the strain is in the form of plastic strain, a macroscopic 
measure of the plastic, permanent and irrecoverable de- 
formations in the system. In this Section we study the 
macroscopic quantities of plasticity such as mean and 
equivalent plastic strains as well as the plastic work, and 
in the next Section in more detail what are the actual 
plastic deformations visible in the void. The concomitant 
dislocations related to void's shape and volume changes 
are studied elsewhere ^L224i 

In deriving the plastic strain here it is assumed that 
the tetragonal symmetry is approximately preserved and 
thus the off-diagonal terms of the stress tensor are neg- 
ligible. Following the literature we separate the strain 
increment de = e dt into elastic and plastic parts^ Thus 
by definition the plastic strain increment becomes 



e%dt 



£af3 dt, 



(9) 
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TABLE I: The onset of plasticity associated with void growth, as indicated by 3 different criteria: deviation from elastic 
behavior in the mean stress, von Mises stress and stress triaxiality. Their threshold values, together with the corresponding 
strain values, are tabulated for uniaxial, biaxial and triaxial expansion at the strain rates e = 10 9 /sec, 5 x 10 8 /sec, 10 s , and 
f0 7 /sec. In particular, the third and fourth columns show the mean stress values a and the corresponding mean engineering 
strain values £ m , respectively, at the critical or yield point at which the mean stress first deviates from the elastic stress-strain 
curve. Analogously, the fifth and sixth columns show the yield point as indicated by von Mises stress a e and the corresponding 
mean engineering strain values e m , respectively. The seventh and eighth columns show the yield point as indicated by the stress 
triaxiality \ an d the corresponding mean engineering strain values e m , respectively. Note the small but significant differences 
in the yield point as indicated by these three different criteria. The error bars of the values are of the order of last reported 
digit. The details of the simulations are in the caption of Fig. 0and the curves, from which the yield data have been calculated, 
are plotted in Figs. 11131 Note that, as expected, von Mises stress is small and erratic in the triaxial case, so those von Mises 
and stress-triaxiality data are not tabulated. 
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where e^°j dt is the total increment of the strain. Below 
we use i instead of i dt since dt can be divided from both 
sides of Eq. JJjJ. The total strain increment is an in- 
put parameter in these strain-controlled simulations. It 
is given by the strain-rate matrix £. The compliance S 
relating the elastic strain increment to the stress incre- 
ment is derived from the stress-strain curves in the elastic 
region by: 

Zapivap) = Sa a p. (10) 

The stress matrix a a p is calculated each time-step using 
Eq. ©. The elastic compliance tensor S(a a p) is retrieved 
from the elastic part of the stress-strain curve of the cases 
without the void as follows. Due to the nonlinearity of a 
stress-strain curve we have not only retrieved the slope 
of it, which would give 3B — Cxi + 2Ci2, where B is 
the bulk modulus, but fitted a fourth order polynomial 
to the strain versus stress curve, whose derivative gives 
us l/ZB~ 1 {a m ). This is done separately for the uniax- 
ial, biaxial, and triaxial no-void cases, and the respective 
curves are used for the cases with the void. It should be 
mentioned, too, that in the derivation of the bulk mod- 
ulus the mean total logarithmic strain is used instead of 
the engineering principal strain used in the plots of this 
Article. Similarly the term C' — \ (Cu — G\%) is derived 
using a fourth order polynomial in the mean strain versus 
von Mises stress curve giving 1/C". Note, that when C 
is derived from the plot using mean strain there are pref- 
actors 1/3 and 2/3 for l/(Cn — C12) in the uniaxial and 



biaxial cases, respectively. Using formulas which relate 
S\x and S12 to Cxx and C12 in cubic crystals^ and the 
correspondence between elastic constants and moduli we 
get for S11 and Si 2 

Using the Su(cr OT ) as S aa {a m ) and Si 2 (cr m ) as S a p(cr m ) 
due to the symmetry and neglecting off-diagonal terms, 
which are small compared to the diagonal ones, we get all 
the necessary terms for S(a m ), and thus E^(u m ) from 
Eq. (I1U|) . Note, that since the a m is used as a parameter 
instead of a a p , the von Mises stress must be mapped with 
the mean stress when finding the corresponding C'. This 
was done again by fitting the von Mises vs. mean stress 
curves with fourth order polynomials. 

Subtracting the elastic strain from the total strain as 
in Eq. © we get the mean plastic strain increment 

= 7T ^ "j £ggi (12) 
a 

which time-integral is plotted in Fig.0]for all the loading 
modes and strain-rates. In these plots one sees that af- 
ter the yield point the mean plastic strain first increases 
roughly exponentially, although the region is too small 
to be definitive, and thenceforth roughly linearly. Note 
that the mean plastic strain is not the equivalent plastic 
strain commonly used in plasticity, which will be defined 
below, but a measure of the porosity. This will be studied 
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in the next Section, where the mean plastic strain will be 
compared with the growth of the volume of the void. 

We turn now to the quantification of the dislocation 
flow, conventionally computed at the continuum level as 
the second invariant of the deviatoric plastic strain, the 
equivalent plastic strain. Typically in the case of tetrag- 
onal total strain, the equivalent plastic strain rate would 
be calculated as: 




1/2 



(13) 



The equivalent plastic strain is calculated in turn as 



ef(*) 



'(t')dt'. 



(14) 



In practice this formula for the equivalent plastic strain 
is problematic in MD for several reasons. First, the time 
and length scales in MD are much shorter than those 
assumed in continuum formulations of plasticity. The 
time scale is a problem because dislocation flow becomes 
partially reversible at short enough time scales. Ther- 
mal fluctuations cause reversible oscillations of disloca- 
tions and fluctuations in the local elastic strain. To the 
contrary, the integrand in Eq. I|14|) is positive definite, 
as appropriate for plastic deformation that is cumula- 
tive even when reversed. In practice, the application of 
Eq. 114fl in MD gives a result dominated by the fluctu- 
ations for small time increments; in fact, in our attempt 
to apply the formula to the MD deformation every 10 
time steps, the contribution of the fluctuations was 22 
times as large as the applied total mean strain (these 
values are obtained from the biaxial case with strain- 
rate e = 10 8 /sec). The formula must be modified to be 
insensitive to thermal fluctuations. Second, the formula 
for the equivalent plastic strain assumes isotropic plas- 
ticity in the following sense. In isotropic plasticity, the 
plastic flow is driven by the shear stress quantified by the 
von Mises stress. The equivalent plastic strain is conju- 
gate to the von Mises stress, and therefore takes on a 
particular signficance in the theory. Implicit is the as- 
sumption, for example, that slip systems that experience 
the same shear stress will exhibit the same plastic strain. 
This assumption is violated in MD for two reasons. Once 
again, the thermal fluctuations may cause the initiation 
of flow on one glide system before that on a symmetri- 
cally related system. This effect is observed in our MD 
simulations. Typically the symmetry is restored after a 
brief period, but because the plastic strain is cumulative, 
the symmetry breaking fluctuation is never eliminated 
from the plastic strain 1|14[) . Second, a more mundane 
reason the assumption of isotropy fails is that the sin- 
gle crystal systems are anisotropic, both because of the 
specific glide planes involved and because of the elastic 
constants, especially in copper. 

It may be possible to rectify these problems while re- 
taining the basic formulation of the equivalent plastic 



strain, for example through a suitable multi-resolution 
calculation of the integral l|14|) . We have made several 
attempts at a new formulation, but we were not able to 
develop a satisfactory algorithm, providing a meaningful 
measure of the plastic strain on MD time scales based on 
the equivalent plastic strain integral l|14|) . We found that 
we could eliminate the anomalies due to fluctuations or 
the anisotropy, but not both simultaneously in a robust 
manner. 

We have therefore turned to a different quantification 
of the plastic strain. Certainly, the full deviatoric plastic 
strain tensor is a measure of the plastic flow, conjugate 
to the deviator stress. Its rate of increase is given by the 
traceless part of Eq. ©. Typically, the rate would be 
integrated in a cumulative fashion, but we will not do so. 
The nature of our simulations is such that at the contin- 
uum level plastic flow is only expected in one direction, 
so any sign reversal may be attributed to fluctuations. 
Wc then calculate 



-a/3 



(*) 



(15) 



where the plastic strain rate is given by Eq. JJjJ. Wc 
emphasize again that no absolute value is taken, so fluc- 
tuations cancel. 

Then in order to have a scalar quantification of the 
plastic strain, we compute the Ji invariant, normalized 
as the equivalent plastic strain would be: 

1/2 




t (t')dt' 



e$ p (t')dt' 



(16) 

We must stress that this quantity is not equal to the 
equivalent plastic strain commonly used in plasticity, ex- 
cept in the extraordinary case of monotonic isotropic 
plasticity. It is not conjugate to the von Mises stress, 
for example, in our MD simulations. Nevertheless, it is a 
useful qualitative measure of the degree of plasticity, and 
it allows us to compare the plastic response as the system 
is loaded in different ways and we call it here equivalent 
plastic strain for simplicity. 

The evolution of the equivalent plastic strain during 
uniaxial and biaxial expansion is plotted in Fig- EI In the 
triaxial case it is essentially zero, as expected by symme- 
try, and therefore it has not been plotted. In practice, 
the stress during triaxial expansion has only a negligi- 
bly small fluctuating shear component, so the calculated 
elastic shear strains are very small, too. Since neither 
the box strain nor the elastic strain has an appreciable 
shear component, the equivalent plastic strain is found 
to be zero. 

Now, once the tensors for both the stress and the plas- 
tic strain are derived, (actually only the diagonal terms 
of the plastic strain are needed) , the plastic work can be 
calculated, 



W P (t) 



(17) 
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(see Fig. |HJ). It should be compared with the tempera- 
ture from the same simulations, Fig. [7J Note, that in 
these simulations, when the dilational strain is applied, 
the thermostat is turned off and thus the temperature is 
allowed to change. First the system cools in the elastic 
regime due to adiabatic cooling on expansion, but when 
plastic deformations begin, work is done in the system 
resulting in heating. We find that the increase in plastic 
work does not match exactly with the temperature. In 
principle we expect several effects to contribute to this 
difference: the surface energy of the void, the defect for- 
mation energies for dislocations and point defects, fur- 
ther adiabatic cooling, and any error in calculating the 
elastic energy or strain from the stress. Using the best 
data available to bound the contributions from surface, 
defect and adiabatic cooling energies, we find that there 
remains an energy deficit that we attribute to an error in 
the calculated elastic energy. The error comes from the 
use of the average stress despite stress inhomogeneity in 
the system due to the void: in the plastic work the prod- 
uct of plastic strain and stress is calculated with averaged 
quantities, while in the temperature the product is cal- 
culated at level of each atom and averaged afterward. 



V. VOID EVOLUTION 

A. Growth of the Void 

We now consider the volume and shape evolution of the 
void. During the MD simulations undergoing expansion, 
the surface of the void is determined by finding individual 
atoms that belong to the surface. This is done by creating 
a fine two-dimensional mesh, in which each mesh point 
corresponds to spherical angular coordinates ((/>, 9) . An 
atom is found to represent the surface at each point of 
the mesh, with some atoms representing multiple mesh 
points. In particular, taking the origin to be the center of 
the void, within the solid angle associated with each mesh 
point, the atom that is closest to the origin is defined 
to be the surface atom at that mesh point. There are, 
however, some uncertainties related to this method. If 
the mesh is too dense with its size diverging it can capture 
almost all the atoms in the system. On the other hand 
if it is too sparse, it may neglect some surface atoms, 
especially when void is anisotropic, non-spherical, and 
has some sharp edges in it. Therefore we introduced a 
width to each of the atoms by drawing a circle around it 
that implies a width (d(j>, d9) to the angles, so that one 
atom can occupy several mesh points in a fine mesh. We 
have typically 75-100 points for each angular coordinate, 
giving a total of 5625-10000 mesh points. In the surface of 
the void there are typically few thousand atoms. Besides 
introducing the width to the atoms, we also select atoms 
based on their radial distance: if an atom has much larger 
radial distance r compared to its neighbors it is neglected 
in order not to capture atoms that do not belong to the 
surface. 



Once the surface atoms are identified, the surface is 
tessellated using a generalization of the Delaunay trian- 
gulation method. 49 The Delaunay triangulation is an op- 
timal triangulation of a collection of points-in our case 
atoms-on the plane. It is optimal roughly in the sense 
that the aspect ratio of the triangles is as near to unity 
as possible; more precisely, the Delaunay theorem guar- 
antees that there is a unique (up to degeneracy) trian- 
gulation such that if any triangle in the triangulation is 
circumscribed by a circle, none of the other points will 
be in the interior of the circle. The Delaunay theorem, 
as formulated, does not apply to points on a curved sur- 
face. In fact, there appears to be a topological obstruc- 
tion to the existence of a unique, optimal triangulation 
on a closed surface when the Euler character is non-zero. 
Nevertheless, it is possible to extend the Delaunay trian- 
gulation algorithm to achieve a locally optimal triangu- 
lation almost everywhere. The approach we have taken 
is to project the points patchwise onto flat surfaces. In 
particular, stereographic projections are used to project 
the upper and lower hemispheres separately onto planes. 
Cylindrical coordinates are used to project the equatorial 
region to a cylinder. The Delaunay algorithm is used to 
triangulate each of these projected regions. The patches 
overlap at latitudes of ±45° where the projections are not 
too distorted. The three patches are sewn together using 
a simple advancing front triangulation at the boundaries. 

Using this triangulation, the volume of the void can be 
calculated precisely by summing up the volumes of the 
tetrahedra with one apex at the center of the void and 
the opposing face on the void surface. As we shall see 
below, the void shape evolves to be far from spherical. 
Therefore the approximation of the void surface by tri- 
angles captures the shape better than just assuming it to 
be spherical and using only the solid angles and radial 
distances of a sphere, when calculating the volume of the 
void. An advantage of this method is that if some atom, 
which should be taken into account, is missing from the 
surface of the void, its position is filled with the triangles 
created by its neighboring atoms, and thus the "hole" is 
well approximated by its neighbors. 

In Fig. [8] the porosity or the ratio between the volume 
of the void and the total volume of the system, 

/ = V void /V, (18) 

is plotted for a fraction of the simulations of the strain- 
rates i = 10 10 /scc, 10 9 /sec, 5 x 10 8 /sec, and 10 8 /sec. It 
should be noted that in order to get information about 
the positions of the atoms for the strain values of the in- 
terest (i.e., close to and after the onset of plasticity) the 
calculations were restarted from already expanded sys- 
tem. After the restart the expansion was applied with 
the same strain-rate as earlier but now to the already 
expanded system, thus the strain-rate was increased by 
a few percent compared to the original [since 7io in Eq. 
@ was the restart value]. Therefore these simulations 
are not precisely comparable with the ones plotted in 
Figs. 17171 where the continuum quantities are shown. In 
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these plots, as for the mean plastic strain, in most cases 
first the void grows exponentially and then (if the cal- 
culation has been carried out that long) it changes to 
a linear growth, see especially Fig. Efb) and the strain- 
rates 10 9 /sec and 5 x 10 8 /sec there. The cross-over to 
the linear growth happens at the same point as the rate 
of decrease of the mean stress slows. Although it is be- 
yond the scope of this Article to go into the analysis, the 
reduction in the growth rate coincides with the point at 
which the dislocation density along the shortest distance 
between the void and its periodic image (at the apices of 
the faceted void, cf. Section lYB|) reaches saturation, and 
the nature of the dislocation activity changes dramati- 
cally. This can be interpreted as a finite size effect as 
the void approaches the boundary of the simulation box 
or as a start of the coalescence process of the void with 
its own periodic image. The void-void interactions and 
the coalescence process of two voids in a less restrictive 
geometry will be presented elsewhere^ The shapes of 
the porosity curves as a whole can be compared with the 
mean plastic strain plotted in Fig. 21 Although the vol- 
ume of the void is not calculated throughout the whole 
simulation, one sees easily, that there is correspondence 
between mean plastic strain and the volume of the void. 
There is of course an offset at the elastic part of the 
simulations, since mean plastic strain equals zero then, 
but the initial volume of the void is finite. The corre- 
spondence will be revisited and studied more carefully in 
Section IV1I However, it can be concluded already here 
that the macroscopic quantity mean plastic strain cap- 
tures the microscopic behavior of the void growth very 
well. Effects such as the excess volume associated with 
defects are negligible. This also means that the matrix 
material is plastically incompressible, the dilation comes 
from the void growth, and thus it is consistent with the 
Gurson type of continuum models. 16 



B. Shape Evolution of the Void 

Let us now look at the shape evolution of the void in 
more detail. In Fig. [§| snapshots of the void are shown 
from uniaxial expansion at the strain-rate e = 10 8 /sec. 
There are several interesting aspects in the snapshots. 
In the first two snapshots at strains e = 5.05% and 
5.26% when the system still behaves elastically, the void 
is expanded in the x-direction, which is the direction of 
the strain. However, after that the void makes a rapid 
shape change and becomes more extended in the trans- 
verse y and z-directions, i.e., the strain-free directions, 
than x-direction. This prolate-to-oblate transition may 
be counterintuitive, but the behavior has been observed 
previously in continuum calculationsj^&SiiS and it has 
been related to the appearance of shallow dimples in the 
fractography studies of ductile fracture surfaces in low 
triaxiality conditions. See also studies of non-spherical 
voidsi^^ For example, Budiansky et al., Ref. l5ll inves- 
tigated void shape change in a non-linear viscous plastic 



flow model. They explained the oblate growth of voids 
under uniaxial loading as due to a non-linear amplifica- 
tion of the shear stress on the surface of the void, with 
the maximal void growth rate at the locations of maxi- 
mal von Mises stress: the equator. Their analysis does 
not apply directly to our simulations since they neglect 
elasticity, and the non-linear viscous solid model they 
have used is not expected to be a precise description of 
the plastic flow early in our MD simulation when the 
prolate-to-oblate transition takes place. Furthermore, it 
is not clear from our simulations what value should be as- 
signed to the strain-rate exponent that controls the non- 
linearity in the model of Budiansky et al, although a 
large value is reasonable. Despite these differences, the 
same localization of plastic flow to the equator of the 
void and the transition to an oblate shape does appear 
in both our simulations and the viscous solid model of 
Budiansky et al. Following some additional expansion 
beyond the transition, the void begins to become faceted, 
as visible in the last three snapshots. It should be men- 
tioned that the anisotropy visible in this uniaxial case 
is less pronounced in the biaxial case. The cases with 
the hydrostatic loading are the most isotropic and the 
octahedral shape, somewhat visible in Fig.^Jf), becomes 
more pronounced^ The octahedral shape has been seen 
in spallation experiments in the FCC metal aluminum, 55 
and also in experiments on the equilibrium void shape of 
another material, silicon, too, where it has been used to 
calculate anisotropic surface energies through an inverse 
Wulff construction^ 6 - In void growth associated with dy- 
namic fracture in copper, several effects contribute to the 
faceting: the low surface energy and high ad-atom energy 
of the {111} surfaces common to FCC metals, the high 
anisotropy of the copper elastic constants (A=3.21), and 
the {111} dislocation glide systems. These effects are 
analyzed in detail elsewhere4£ 

In order to characterize the shape change of the void 
not only qualitatively and visually, but also quantita- 
tively, multipole moments of the void shape are calcu- 
lated. To the best of our knowledge, this is the first time 
that multipole moments have been used to characterize 
surface shape. They are a powerful way to quantify the 
evolution of the complex surface containing thousands of 
atoms, and they are suitable for use in continuum mod- 
els and experimental void characterization as well. Using 
spherical harmonics, 

Y lm (f)=Y lm {e(r),4>(f)}, (19) 

expressed as polynomials of Cartesian coordinates, in 
contrast to more commonly used trigonometric forms^i 
we are able to define different multipole moments of the 
void based on its surface atoms: 

Qim = £ fY lm {6,cj ) )r 2 (6,(b)dn, (20) 

where the mean square radius f 2 = j- J r 2 (#,</>) dil. This 
is in contrast to the volume integral more commonly used 
when calculating multipole moments. The axial index of 
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the moment ranges m = and for each m except 

m = 0, Qi m has both real and imaginary parts. Here 
we concentrate on I = 1,2,3, and 4. Only the positive 
moments of m are calculated, since the negative ones are 
related by 



Qi,- 



(21) 



In all 24 different terms are calculated. The polynomial 
forms used here of the moments are listed in the Ap- 
pendix. 

The moments Qi m are not rotationally invariant, but 
depend on the way the coordinates x, y, and z are cho- 
sen. The set of (21 + 1) moments at fixed I form an 
irreducible representation of the SO (3) rotation group, 
and are mixed by rotations according to the usual trans- 
formation rules. They may be combined into a single 
rotationally invariant combination for each I according 
to 



Qi = 



21 



1 - 



m— — l 



1/2 



(22) 



see e.g. Ref. 59. Their use drastically reduces the amount 
of data to be shown. Only the positive m's are needed 
for Qi due to the square of the norm of Qi m and Eq. JUJ. 

Technically the calculation of the multipolc moments 
has been done using the information about the shape of 
the void obtained from the surface triangulation proce- 
dure explained earlier. As in the calculation of the void 
volume, some refinements have been introduced to re- 
duce the uncertainty in the values of the moments that 
arises from single atoms near the threshold for inclusion 
as surface atoms. These borderline atoms can appear in- 
termittently on the void surface during the growth, and 
the tessellation is used to minimize their effect on the 
moments. In calculating the volume of the void the tri- 
angulation gave one face of the tetrahedra that acted as 
small volume elements for the total volume. Here the tri- 
angulation is used to weight the atoms by the amount of 
solid angle associated with each surface atom. In particu- 
lar, each triangle in the tessalation contributes one third 
of its projected solid angle to each of the three atoms that 
make up its vertices, where the solid angle of a triangle is 
computed using the formula 8fl = A1+A2+A3 — tt, where 



.4 



-Ci+lC i + 2 



— — . and Ci = Xi+i ■ x t+2 for 

V i+lV " 

i = 1,2,3 (mod 3) and where xi is the unit vector in the 
direction of the i th vertex of the triangle.— The weight of 
each atom is the sum of these solid angle contributions. 
This reduces the sensitivity of moments to the atomic 
discretization of the surface, since evanescent atoms that 
occasionally appear and disappear from the fluctuating 
surface only make a small, local change to the value of 
r 2 . It may be of interest to note that in the course of the 
development of these surface multipolc moments, several 
different variations on the definition of the moments were 
tried. The definition presented here (|20|l produced sub- 
stantially less noise (up to a factor of 5 less noise) than the 



other definitions we tried, even though they all showed 
the same trends in the shape evolution. Using these 
weights for the atoms and after normalizing each atom's 
(x, y, z) coordinate by its distance r = (x 2 + y 2 + z 2 ) 1 / 2 
from the center of the void all the terms in Eqs. (|A1|) - 
(|A4(I , and Eq. \\12l , are calculated. The center of the void 
is defined to be the point where the three components of 
Qim, as given by Eq. IjAlfl . are zero. 

Due to space limitations, only a fraction of the mul- 
tipole moment data is shown here. In Fig. HOf a) the 
quadrupole moments Q2m for all the positive m of the 
void are shown in the uniaxial case for the strain-rate 
i = 10 8 /sec. This is the same simulation as the snapshots 
in Fig. Indeed, the quadrupole moments are able to 
represent numerically the shape changes one sees in the 
snapshots. In Figs.^a) and (b) at strains e = 5.05% and 
5.26% the void is elongated to the direction of the load, 
which is visible as Q20 > 0. Between strains e = 5.47% 
[Fig- Etc)] and 5.68% [Fig.HJd)] the void is extended more 
transverse to the direction of the strain, thus Q20 < 0, 
and later it starts to becomes more of octahedral shape 
and the absolute value of Q2 saturates. 

Figures HUT bWd) show the rotational invariant multi- 
pole moments Qi, Eq. (|22|l . in cases with uniaxial, biax- 
ial, and triaxial loading, respectively. Each of the cases 
has strain-rate e = 10 8 /sec. In the plots it is clear that 
the behavior that the quadrupole moment has first a non- 
zero value and then makes a rapid dip but returns back 
to a non-zero value due to the transverse elongation is the 
strongest in the uniaxial case. On the other hand the oc- 
tahedral shape measured by the Q4 is more pronounced 
in the biaxial and triaxial cases than in the uniaxial case 
as explained qualitatively earlier. Hence we find that the 
multipole moments introduce a good method to measure 
the shape changes of the void. The non-zero values for 
Q3, as well as Q2 in other cases than uniaxial, indicate 
that the void is not (cubically) symmetric in these simu- 
lations. It should be mentioned that these first four mo- 
ments are enough to characterize the shapes of the void 
and the higher moments contain little relevant informa- 
tion. This was checked by creating a three dimensional 
surface based on the moment values and drawing it in the 
same figure with the actual positions of the surface atoms 
using a standard commercial visualization program. The 
surfaces overlapped very well. 



VI. SUMMARY OF THE UNIAXIAL CASE 

Based on the data shown earlier in this Article for 
the shape and volume changes of the void as well as the 
stress-strain behavior and the stress-triaxiality, it is ev- 
ident that the uniaxial loading raises many interesting 
aspects to be studied in more detail. Therefore we now 
concentrate on the uniaxial case when summarizing how 
the void evolves and how the evolution is related to the 
stress-triaxiality as the system is expanded. By plotting 
most of the measured values together in one figure for 
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the uniaxial simulation at the strain-rate e = 10 8 /sec, 
it is possible to compare the evolution sequence and in- 
vestigate causality, see Fig. UTT a^. For clarity, we have 
chosen not to plot many quantities in the figure, e.g. 
plastic strain, plastic work, and temperature. However, 
their concomitant behaviors are included in the verbal 
explanation below and shown in previous figures. The 
data shown in this figure are from the restarted simula- 
tion (see the explanation near Fig. |SJ, as are the data 
in Figs. 1811171 The mean stress and stress-triaxiality data 
are from that simulation, too, and thus arc different from 
the data shown in Figs.^a) and GJa). In any case, the 
overall behavior stays the same as well as the other de- 
tails as the system size, etc. 

The evolution of the void and the system's stress-strain 
behavior can be divided in three or even four different 
regions. The first region is when the system expands 
elastically. The mean and von Mises stresses increase 
smoothly, nearly linearly, and the stress-triaxiality stays 
nearly constant. Through the elastic region the void vol- 
ume fraction remains nearly constant, too. It is not ex- 
actly constant, since due to the free surface of the void, 
the elastic expansion is a bit greater at the surface of the 
void compared to the total system. Trivially the mean 
and equivalent plastic strains as well as the plastic work 
are equal to zero in the elastic region, and temperature 
decreases in the system. The quadrupole moment has a 
non-zero value because of the elongation in the direction 
of the strain. 

The second region begins at what we call the yield 
point, i.e., the onset of rapid growth of the void facili- 
tated by plastic deformation. Heterogeneous nucleation 
of dislocations at the void surface is the primary mech- 
anism for plasticity in the simulation, and thus it is at 
this point that the measured quantities start to deviate 
from their elastic behavior. The mean stress begins to 
plateau here, but fluctuating somewhat. The early depar- 
ture from elastic behavior prior to the plateau is much 
less pronounced. The change in the void shape begins 
just at the point the mean stress deviates from elastic 
behavior: Q20 goes rapidly from the positive value ac- 
quired during elastic expansion to a negative value, i.e., 
from a prolate shape (elongated in the direction of the 
strain) to an oblate shape (expanded in the transverse 
directions). Q2, on the other hand, drops from a positive 
value, almost reaching zero at the prolate-to-oblate tran- 
sition point (e = 5.45%) and rising even larger value after 
that [in fact, the oblate shape is somewhat more pro- 
nounced than the earlier prolate shape, seen as a larger 
absolute value of Q20 in Fig. HUt a - )]. When Q2 starts to 
drop Qa starts to rise. Then after the prolate-to-oblate 
transition point, Q4 begins to saturate. At a strain of 
e = 5.55%, Q2 is 1.5 times as large as its value at the 
end of the elastic phase (e = 5.25%). Mean plastic strain, 
equivalent plastic strain, plastic work, and temperature 
increase together with porosity. Their increase starts im- 
mediately at the yield point, i.e., when the mean stress 
first deviates from the elastic behavior. A bit later than 



the plastic strain, equivalent plastic strain, plastic work, 
and temperature, the stress triaxiality increases simulta- 
neously with the first substantial drop in the von Mises 
stress. The fact that the increase in stress-triaxiality fol- 
lows later is dependent on how the ratio between mean 
stress and von Mises stress develops, as discussed ear- 
lier. In Fig. Ota) and in the data reported in Table [I] the 
stress-triaxiality started to increase simultaneously with 
the mean and von Mises stresses deviating from the elas- 
tic behavior. The increase of stress-triaxiality is caused 
by von Mises stress plummeting in contrast to nearly con- 
stant mean stress. The drop in von Mises stress follows 
from the flow of dislocations nucleated at the void and 
from the relaxation of the shear stress of the system due 
to the flow. 

The third region is when the void fraction, mean plastic 
strain, equivalent plastic strain, plastic work, and tem- 
perature switch from rapid increase to a linear growth or 
even saturate. Subsequently the increase of the stress- 
triaxiality slows down and plateaus. The value at the 
plateau is related in continuum models to the ratio of 
the mean stress threshold for void growth to the flow 
stress. At the plateau von Mises stress saturates at a 
small value, close to the tensile strength, and the shape 
of the void starts to become more of octahedral shape 
although having a non-zero quadrupole moment, too. 
Hence at the second and third regions the mean stress 
is nearly constant, but von Mises stress and the stress- 
triaxiality changes. 

A conclusion might be that once the threshold for void 
growth is reached, the population of dislocations rises 
sufficiently to relax the shear stress quite effectively and 
it drops to a low value (the flow stress); the mean stress, 
on the other hand, plateaus since it is relaxed by void 
growth and requires that the stress at the void surface 
be sufficiently high to continue to nucleate dislocations. 
The fourth region is the failure, when the system breaks, 
and it is not studied here. 

In order to see if the rapid changes studied above are 
due to the smallness of the size of the void we have done 
one additional simulation with a system in which the 
initial void radius is half that in the other simulations; 
otherwise, the system size is the same, see Fig- lllf b 1 ). For 
the small void simulation, uniaxial strain at a strain-rate 
i = 10 8 /sec has been used. Here the difference is that the 
quadrupole moment suffers stronger fluctuation due to 
the smallness of the void, where each of the surface atom 
contribute more to its value and therefore is even more 
sensitive to the selection criterion of surface atoms, and 
also the shape changes are harder to determine based on 
the Q20 behavior. The other main difference is that the 
growth of the void is linear all the time. Also the changes 
in stress-triaxiality in fact advances the porosity when 
saturating, and the mean stress does not fluctuate but 
drops more rapidly (this can be compared with the case 
without the void, where the mean stress drops abruptly) . 

We have also compared the mean plastic strain and 
the void volume fraction calculations. In continuum solid 
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mechanics, it is assumed that solid materials are plasti- 
cally incompressible. Any local dilation, as indicated by 
a change in the mean strain, is attributed either to elastic 
dilation or to a change in the porosity of the material, 
where the porosity is equated to the void volume frac- 
tion. The porosity / and the mean plastic strain are 
then related according to the equation^ 

/ = 3(l-/)«£ (23) 

where the dots denote time derivatives. Integration with 
respect to time, porosity from /o to /, and mean plastic 
strain from zero to e^, gives 

/(4) = l + (/o-l)exp(-3 £ Q. (24) 

where fo is the initial porosity. It is interesting to check 
whether this relationship holds for the MD simulation, 
where other effects such as excess volume associated with 
dislocation cores, vacancies or other defects could require 
corrections. In comparing the porosity inferred from the 
mean plastic strain and that calculated directly as the 
void fraction, the agreement is very good. The trends 
are in excellent agreement, but there is a small discrep- 
ancy between the curves, so that the porosity from the 
mean plastic strain is overestimated. We believe that the 
discrepancy arises because of the void surface. In calcu- 
lating the void fraction, we have defined the void surface 
to pass through the center of the surface atoms. How- 
ever, the properties of the surface atoms are distinct from 
the bulk atoms. Therefore, there is some ambiguity in 
where the surface should be placed, and a small uniform 
shift 5r of the surface radially into the bulk is enough 
to account for the discrepancy. In Fig. ^] we have plot- 
ted the comparison of the porosity from the mean plastic 
strain, Eq. Q24JI. and from the void fraction, Eq. <|18() . us- 
ing a constant radius increase Sr = 0.58 ag, where ag is 
the lattice constant for the void volume calculation. The 
correction for the void size, Sr, is a fit parameter and it 
varies for different strain-rates and slightly for different 
loading modes, but is always positive and of the order 
of the lattice constant, do- It should be noted, too, that 
by Taylor expanding Eq. I|24|) and discarding the higher 
order terms, it becomes 

/(4) = /o + 3(l-/ )C (25) 

indicating a linear correspondence between / and e^ n as 
long as the void fraction is small. 

VII. CONCLUSIONS AND DISCUSSIONS 

In this Article void growth in copper has been stud- 
ied in a high range of strain-rates at the atomistic level. 
The model has been designed to simulate the growth 
of a void nucleating from a very weakly bound inclu- 
sion during strain-controlled dynamic fracture. In or- 
der to see the effect of various modes of expansion and 



the related stress-triaxiality, three different modes have 
been applied, namely uniaxial, biaxial, and triaxial. The 
molecular dynamics method developed here has been 
shown to be efficient enough to explore the different load- 
ing conditions and strain-rates varying over four orders 
of magnitude. A uniform expansive loading of a sys- 
tem with periodic boundary conditions has been imple- 
mented using a well defined scaling matrix method. For 
the longest calculations, the MD method was parallelized 
successfully. The macroscopic observables mean stress, 
von Mises stress, stress-triaxiality, mean plastic strain, 
equivalent plastic strain, plastic work, and temperature 
have been calculated and compared with the microscopic 
quantities measured at the atomistic-level, such as the 
volume of the void and its shape change. A method to 
describe the shape changes in the void is introduced and 
employed, namely calculations of the multipole moments 
of the void based on spherical harmonics in polynomial, 
not trigonometric, form. When calculating the volume of 
the void with an unknown shape or defining solid surface 
for the multipole moment calculation a useful method, 
namely optimal triangular tessellation, has been intro- 
duced. This method has been extended from the usual 
planar case to non-planar objects such as the surface of 
the void. 

When the different measured quantities are compared 
with each other during an MD simulation in uniaxial ex- 
pansion, it is found that at early stages of plasticity von 
Mises stress, and thus also stress-triaxiality, plays a more 
significant role to the void growth and its shape change 
than expected. On the other hand, most of the macro- 
scopic plastic quantities as mean and equivalent plastic 
strain as well as plastic work and temperature, seem to 
be more dependent on the simultaneous saturation of the 
mean stress. These calculations show a counter-intuitive 
behavior, observed previously in continuum void growth 
modelingj22i£ii££ that a prolate-to-oblate transition oc- 
curs. When the system starts to yield, the expansion of 
the void switches from its original elastic extension in the 
direction of the load to transverse plastic expansion. 

The yield stress values for the lowest strain-rates 
10 7 /sec are in reasonable agreement with the experimen- 
tally measured spall strength^ The fact that mean plas- 
tic strain can be mapped to the growth of the void is 
consistent with continuum models^ 

This study leaves many open questions. For instance 
related to the void growth are the dislocations, which 
occur when the system yields. Since the FCC crystal 
studied here is perfect apart from the void, the disloca- 
tions form from void's surface. They are also responsi- 
ble for its growth by carrying material away. Thus the 
characterization of plasticity surrounding a growing void 
at the level of dislocations should be investigated, too, 
especially with respect to the stress-triaxiality. These in- 
vestigations are underway4£ Their results are beyond the 
scope of this Article, other than to mention that the iden- 
tification of the yield point in this Article does indeed cor- 
respond to the point of initial nucleation of dislocations. 
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Another topic that is beyond the scope of this Article 
and needs further investigation, but is closely related to 
the studies here, is the quantitative connection between 
the shear stress, and thus the mode of the loading, and 
the onset of the void growth and the resulting change in 
the stress state. Other areas where this study can be ex- 
tended are different materials including different lattice 
structures such as body-centered cubic (BCC) lattices; 
in the uniaxial case other orientations of the lattice as 
(110) and (111); continuously changing stress-triaxiality 
in order to create the full yield surface to the stress space 
similarly as in Gurson type of continuum studies^ to in- 
clude grain boundaries, defects, pre-existing dislocations, 
several voids, etc. 
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APPENDIX A: MULTIPOLE MOMENTS 



Terms with / = 2 are the quadrupole moments and 
they get non-zero values if there is ellipsoidal shape in 
the object. The quadrupole moments are as follows: 

Q20 = l^^j3z 2 ~r 2 dfi, 
ReQ2i = -^ v /|f^ Jxzdfl, 
Im Q 21 = -±^£^JyzcM, (A2) 

Im Q 2 2 = ky/Ew Jxydn. 

Terms with I = 3 are the octupole moments and they 
get non-zero values for tetrahedron shapes: 

Re Q 31 = -i^/f £ j>(5z 2 - r 2 ) dfi, 

Im Q 31 = -lJ^^Jly(5z 2 -r 2 )<m, 

ReQ 32 = ±^^J±z(x 2 -y 2 )dn, ( A 3) 

Im Q 32 = l^^jl xyzd n, 

Re Q33 - -kJfw IH x3 ~ 3xy 2 )dn, 

Im Q 33 = |^/f ^/i^-S^dfi. 



The 24 different polynomial terms of the multipole mo- 
ments used in this study are listed below m& Here the con- 
ventional notation is used, so that the principal axis of 
the coordinates is the z-direction. When these formu- 
las are used in interpreting the void shape evolution, the 
principal axis is the direction with uniaxial loading, i.e., 
the x-coordinate in the Article. Similarly the load free 
directions y and z correspond to x and y below. 

The polynomial terms when I = 1 are the dipole mo- 
ments and they capture if the object is offset. The dipole 
moments are: 

Re Qn = jrxdVL, ^ 

Im Q n = -ly^LJ^ Jrydn, 

where f 2 = ± J r 2 (d, </>)dfi. 



And finally the terms with I — 4 are listed. They 
are the hexadecapole moments and capture octahedron 
shapes: 

Q40 = J>(3r 4 - 30r 2 z 2 + 35z 4 ) dfi, 

Re Q41 = J ^xz(7z 2 - 3r 2 ) dfi, 

Im Q41 = -| y/% £ / ^yz(7z 2 - 3r 2 ) dfi, 

Re Q42 = l/i^JM* 2 ~ y T K 7z2 - rT ) dSl , 

Im Q 42 - f±xy(7z 2 - r 2 ) dfi, (A4) 

Re Q43 = -f ^/fw I M x3 - 3xy 2 )zdn, 

im q 43 = -| yjf ^JM 3x2 y - y 3 ) zdn - 

Re Q44 = I^(x 4 - 6x 2 y 2 + y 4 ) dfi, 

im Q44 = f^f w J^xyix 2 - y 2 ) dn - 
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FIG. 3: The stress-triaxiality \ J7| versus engineering strain 
e from the same simulations as in Fig. Q See the caption of 
Fig-0for the details, (a) Uniaxial and (b) biaxial expansion. 
In triaxial expansion stress-triaxiality is diverging and not 
defined. 



FIG. 4: Mean plastic strain e^, calculated using Eq. 1121 . 
versus engineering strain e from the same simulations as in 
Fig. See the caption of Fig. for the details, (a) Uniaxial, 
(b) biaxial and (c) triaxial expansion. 



FIG. 5: Equivalent plastic strain ef , calculated as Eq. I16L 
versus engineering strain e from the same simulations as in 
Fig. Q See the caption of Fig. Q for the details, (a) Uniaxial 
and (b) biaxial expansion. 



FIG. 6: Plastic work Wp, calculated from Eq. 1171 . versus 
engineering strain e from the same simulations as in Fig. Q 
See the caption of Fig. for the details, (a) Uniaxial, (b) 
biaxial, and (c) triaxial expansion. 



FIG. 7: Temperature T, versus engineering strain e from the 
same simulations as in Fig. Compare with the plastic work 
plotted in Fig.|H] See the caption of Fig.Qfor the details, (a) 
Uniaxial, (b) biaxial, and (c) triaxial expansion. 



FIG. 8: Void volume fraction versus engineering strain e. The 
evolution of the ratio of the void volume to the total box 
volume is plotted for strain-rates e = 10 10 /sec, 10 9 /sec, 5 x 
10 8 /sec, and 10 8 /sec. See the caption of Fig.Qfor additional 
details, (a) Uniaxial, (b) biaxial, and (c) triaxial expansion. 



FIG. 9: Snapshots of the atoms comprising the surface of the 
void during uniaxial expansion with e x — e, e y — e z — 0. 
The simulation box is oriented along (100) directions, so that 
the z-axis is out of the paper. The strain-rate is e — 10 8 /sec. 
See the caption of Fig. 0^a) for additional details. The panels 
show snapshots at different strains: (a) e = 5.05% (b) e = 
5.26% (c) e = 5.47% (d) e = 5.68% (e) e = 5.89% (f) e = 
6.10% 



FIG. 10: Multipole moments of the void surface calculated 
using Eq. 120H . (a) Quadrupole moment Qim with m — 0, 1, 2 
for uniaxial expansion at e = 10 8 /sec. (b)-(d) The moments 
Qi for I = 1, 2,3,4 in (b) uniaxial, (c) biaxial, and (d) triaxial 
cases. They are calculated using Eq. 122H and the strain- 
rate i = 10 8 /sec. See the caption of Fig. 0for details of the 
simulations. 
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FIG. 11: (a) The mean stress a m (thick solid line), stress- 
triaxiality % (dotted line), volume fraction / of the void 
(dashed line), and the quadrupole moment Q20 (thin solid 
line) from the simulation with uniaxial expansion at e = 
10 s /sec. See the captions of Figs. □ El [HI and [TO] for the 
details, (b) As a comparison the same measures as in (a), but 
now for the case having an initial void radius of 1.1 nm and 
863 543 atoms in the system undergoing uniaxial expansion 
at the same e = 10 s /sec. 



FIG. 12: Porosity / calculated from the actual void fraction 
as in Eq. 1181 with 5r — 0.58 ao (see the text for details of Sr) 
and from the mean plastic strain as in Eq. 1241 from the 
simulation with uniaxial expansion at i = 10 8 /sec. See the 
captions of Figs. HI and llll for the simulation details. 
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